1 PCAngsd Selection Scan

  • Genome-wide Fst scan between years of PWS pop using PCAngsd

1.1 Run PCAngsd selection on Pairwise comparison

  • Can be run for multiple populations
python /home/jamcgirr/apps/pcangsd/pcangsd.py -beagle /home/ktist/ph/data/new_vcf/MD7000/beagle/subpop/PWS91.PWS96_maf05_BEAGLE.PL.gz -o /home/ktist/ph/data/angsd/selection/PWS91.PWS96_selection -selection -sites_save 

# scripts to run at Fram: pcangsd_selection.sh 
# run pcansgd_selectionPWS.sh for aggregated PWS pop over years

1.2 Results

cols<-c("#0072b2","#cc79a7","#009e73","#d55e00","#56b4e9","#e69f00","#f0e442")

pop_info<-read.csv("../Data/Sample_metadata_892pops.csv")
pop_info<-pop_info[,c("Sample","pop","Year.Collected")]
colnames(pop_info)[3]<-"year"

#Based on the tutorial http://www.popgen.dk/software/index.php/PCAngsdTutorial

## function for QQplot
qqchi<-function(x,...){
    lambda<-round(median(x)/qchisq(0.5,1),2)
    qqplot(qchisq((1:length(x)-0.5)/(length(x)),1),x,ylab="Observed",xlab="Expected",...);abline(0,1,col=2,lwd=2)
    legend("topleft",paste("lambda=",lambda))
}

########  PWS pairwise run ###
pwapops<-c("PWS91" ,"PWS96","PWS07", "PWS17")
comb<-combn(pwapops, 2)
comb<-t(comb)
comb1<-comb[c(1,4,6,3),]

for (i in 1:nrow(comb1)){
    pop1<-comb1[i,1]
    pop2<-comb1[i,2]
    
    #s<-npyLoad(paste0("Data/PCAangsd/selection/",pop1,".",pop2,"_selection.selection.npy"))
    s<-npyLoad(paste0("../Data/PCAangsd/selection/pruned_",pop1,".",pop2,"_selection.selection.npy"))
    
    #how many PC axes were evaluated?
    nc<-ncol(s)
    
    ## make QQ plot to QC the test statistics
    #qqchi(s)
    
    # convert test statistic to p-value
    if (nc==1) pval<-1-pchisq(s,1)
    if (nc>1) {
        p$pval1<-pval[,1]
        p$pval2<-pval[,2]
        p$loc<-1:nrow(p)
        p$pval1.log<--log10(p$pval1)
        p$pval2.log<--log10(p$pval2)
        }
    
    #read the position info    
    p<-read.table(paste0("../Data/PCAangsd/selection/pruned_",pop1,".",pop2,"_selection.sites"),colC=c("factor","integer"),sep=":")
    names(p)<-c("chr","pos")
    
    ## make manhatten plot
    p$pval<-pval
    p$loc<-1:nrow(p)
    p$pval.log<--log10(p$pval)
    
    write.csv(p, paste0("../Output/PCA/selection/Fst_pval_Pruned", pop1,".",pop2,".csv"))
    write.csv(p[p$pval.log>4,], paste0("../Output/PCA/selection/Selected_highP_sites_Pruned", pop1,".",pop2,".csv"))
    
    ch<-unique(p$chr)
    
    #count the number of sites per chromosomes
    poss<-data.frame(chr=ch)
    k=1
    for (i in 1:nrow(poss)){
        df<-p[p$chr==poss$chr[i],]
        poss$start[i]<-k
        poss$end[i]<-k+nrow(df)-1
        k=k+nrow(df)
    }
    
    poss$x<-poss$start+(poss$end-poss$start)/2
    
    #color vectors
    colors<-rep(c("steelblue","lightblue"), times=13)
    p$chr<-factor(p$chr, levels=poss$chr)
    ggplot(data=p, aes(x=loc, y=pval.log, color=chr))+
        geom_point(size=0.15)+
        scale_color_manual(values=colors, guide='none')+
        scale_x_continuous(name="Chromosome", breaks=poss$x, labels=gsub("chr",'',poss$chr))+
        theme_classic()+ylab("-log10(p-value)")+
        ggtitle(paste0(pop1," vs. ",pop2))
    ggsave(paste0("../Output/PCA/selection/pcangsd_selection_",pop1,"_",pop2,"_plot.png"), width = 6, height = 3, dpi = 300)    
}


1.3 PCAnsgd for all years of PWS

## PWS together
s<-npyLoad("../Data/PCAangsd/selection/PWS_selection.selection.npy")

## make QQ plot to QC the test statistics
#qqchi(s)
# convert test statistic to p-value
pval<-1-pchisq(s,1)

## read positions 
p<-read.table("../Data/PCAangsd/selection/PWS_selection.sites",colC=c("factor","integer"),sep=":")
names(p)<-c("chr","pos")
p$pval1<-pval[,1]
p$loc<-1:nrow(p)
p$pval1.log<--log10(p$pval1)

#count the number of sites per chromosomes
ch<-unique(p$chr)
poss<-data.frame(chr=ch)
k=1
for (i in 1:nrow(poss)){
    df<-p[p$chr==poss$chr[i],]
    poss$start[i]<-k
    poss$end[i]<-k+nrow(df)-1
    k=k+nrow(df)
}

poss$x<-poss$start+(poss$end-poss$start)/2

colors<-rep(c("steelblue","lightblue"), times=13)
p$chr<-factor(p$chr, levels=poss$chr)

ggplot(data=p, aes(x=loc, y=pval1.log, color=chr))+
    geom_point(size=0.1)+
    scale_color_manual(values=colors, guide='none')+
    scale_x_continuous(name="Chromosome position", breaks=poss$x, labels=gsub("chr",'',poss$chr))+
    theme_classic()+ylab("-log10(p-value)")+
    ggtitle("PWS")
ggsave("../Output/PCA/selection/PWS_pcangsd_selection.png", width = 12, height = 6, dpi=300)    

1.4 High P-value loci from each pop pair

#for each population pair
pwapops<-c("PWS91" ,"PWS96","PWS07", "PWS17")
comb<-combn(pwapops, 2)
comb<-t(comb)
comb1<-comb[c(1,4,6,3),]

#create windows to be assigned for each site
chr <- read.table('../Data/new_vcf/chr_sizes.bed')
chr<-chr[,-2]
colnames(chr)<-c("chr", "len")
#chrmax$start<-c(0,cumsum(chrmax$len)[1:(nrow(chrmax)-1)])
#setkey(chrmax, chr)

winsz <- 50000 # window size
# use seq to find the start points of each window
windows<-data.frame()
for (i in 1:nrow(chr)){
    ch<-chr$len[chr$chr==paste0("chr",i)]
    window_start <- seq(from = 1, to = ch-50000, by = winsz)
    window_stop <- window_start + (winsz-1)
    win <- data.frame(start = window_start, stop = window_stop, 
                      mid = window_start + ((window_stop-window_start)/2-0.5))
    win$chr<-paste0("chr",i)
    win$ch<-i
    windows<-rbind(windows, win)
}
setDT(windows)
fstdat<-list()
Fst<-data.frame()
for (i in 1: nrow(comb1)){
    pop1<-comb1[i,1]
    pop2<-comb1[i,2]
    df<-fread(paste0("../Output/PCA/selection/Selected_highP_sites_Pruned", pop1,".",pop2,".csv"))
    df<-df[,-1]
    #assign the window info
    df<-windows[df, .(chr, pos, pval, pval.log, start, stop, mid), on=.(chr=chr, start <=pos, stop>=pos)]
    #create a unique id for the window
    df$win_id<-paste0(df$chr,"_",df$mid)
    df$pops<-paste0(pop1,".",pop2)
    fstdat[[i]]<-df
    names(fstdat)[i]<-paste0(pop1,".",pop2)
    Fst<-rbind(Fst, df)
}

#Find overlapping positions among pairs?
lapply(fstdat, "[[", "win_id")

Reduce(intersect, list(fstdat[[1]]$win_id , fstdat[[2]]$win_id))
Reduce(intersect, list(fstdat[[2]]$win_id , fstdat[[3]]$win_id))
Reduce(intersect, list(fstdat[[3]]$win_id , fstdat[[4]]$win_id))

Fst$chr<-factor(Fst$chr, levels=paste0("chr",1:26))
ggplot(Fst, aes(x=pos, y=pval.log, color=pops))+
    geom_point()+
    facet_wrap(~chr)+
   theme_bw()+xlab("Genome position")+
    ylab("-log10(P-val)")+theme(legend.title = element_blank(), axis.text.x = element_blank())
ggsave("../Output/Fst/pcangsd_scan_highPsites.png", width=9, height=6, dpi=300)



LS0tCnRpdGxlOiAiU2VsZWN0aW9uIHNjYW4gKEZzdCBvdXRsaWVycykgYmV0d2VlbiB5ZWFycyIKb3V0cHV0OgogIGh0bWxfbm90ZWJvb2s6CiAgICAgIHRvYzogdHJ1ZSAKICAgICAgdG9jX2Zsb2F0OiB0cnVlCiAgICAgIG51bWJlcl9zZWN0aW9uczogdHJ1ZQogICAgICB0aGVtZTogbHVtZW4KICAgICAgaGlnaGxpZ2h0OiB0YW5nbwogICAgICBjb2RlX2ZvbGRpbmc6IGhpZGUKICAgICAgZGZfcHJpbnQ6IHBhZ2VkCi0tLQpgYGB7ciBldmFsPUZBTFNFLCBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFLCBpbmNsdWRlPUZBTFNFfQpzb3VyY2UoIi4uL1JzY3JpcHRzL0Jhc2VTY3JpcHRzLlIiKQpsaWJyYXJ5KHN0cmluZ3IpCmxpYnJhcnkoZ3JpZEV4dHJhKQpsaWJyYXJ5KFJjcHBDTlB5KQoKYGBgCgojIFBDQW5nc2QgU2VsZWN0aW9uIFNjYW4KKiBHZW5vbWUtd2lkZSBGc3Qgc2NhbiBiZXR3ZWVuIHllYXJzIG9mIFBXUyBwb3AgdXNpbmcgUENBbmdzZCAgIAoKIyMgUnVuIFBDQW5nc2Qgc2VsZWN0aW9uIG9uIFBhaXJ3aXNlIGNvbXBhcmlzb24KKiBDYW4gYmUgcnVuIGZvciBtdWx0aXBsZSBwb3B1bGF0aW9ucyAgCgpgYGB7YmFzaH0KcHl0aG9uIC9ob21lL2phbWNnaXJyL2FwcHMvcGNhbmdzZC9wY2FuZ3NkLnB5IC1iZWFnbGUgL2hvbWUva3Rpc3QvcGgvZGF0YS9uZXdfdmNmL01ENzAwMC9iZWFnbGUvc3VicG9wL1BXUzkxLlBXUzk2X21hZjA1X0JFQUdMRS5QTC5neiAtbyAvaG9tZS9rdGlzdC9waC9kYXRhL2FuZ3NkL3NlbGVjdGlvbi9QV1M5MS5QV1M5Nl9zZWxlY3Rpb24gLXNlbGVjdGlvbiAtc2l0ZXNfc2F2ZSAKCiMgc2NyaXB0cyB0byBydW4gYXQgRnJhbTogcGNhbmdzZF9zZWxlY3Rpb24uc2ggCiMgcnVuIHBjYW5zZ2Rfc2VsZWN0aW9uUFdTLnNoIGZvciBhZ2dyZWdhdGVkIFBXUyBwb3Agb3ZlciB5ZWFycwoKYGBgCgojIyBSZXN1bHRzCgpgYGB7ciBldmFsPUZBTFNFLCBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQpjb2xzPC1jKCIjMDA3MmIyIiwiI2NjNzlhNyIsIiMwMDllNzMiLCIjZDU1ZTAwIiwiIzU2YjRlOSIsIiNlNjlmMDAiLCIjZjBlNDQyIikKCnBvcF9pbmZvPC1yZWFkLmNzdigiLi4vRGF0YS9TYW1wbGVfbWV0YWRhdGFfODkycG9wcy5jc3YiKQpwb3BfaW5mbzwtcG9wX2luZm9bLGMoIlNhbXBsZSIsInBvcCIsIlllYXIuQ29sbGVjdGVkIildCmNvbG5hbWVzKHBvcF9pbmZvKVszXTwtInllYXIiCgojQmFzZWQgb24gdGhlIHR1dG9yaWFsIGh0dHA6Ly93d3cucG9wZ2VuLmRrL3NvZnR3YXJlL2luZGV4LnBocC9QQ0FuZ3NkVHV0b3JpYWwKCiMjIGZ1bmN0aW9uIGZvciBRUXBsb3QKcXFjaGk8LWZ1bmN0aW9uKHgsLi4uKXsKICAgIGxhbWJkYTwtcm91bmQobWVkaWFuKHgpL3FjaGlzcSgwLjUsMSksMikKICAgIHFxcGxvdChxY2hpc3EoKDE6bGVuZ3RoKHgpLTAuNSkvKGxlbmd0aCh4KSksMSkseCx5bGFiPSJPYnNlcnZlZCIseGxhYj0iRXhwZWN0ZWQiLC4uLik7YWJsaW5lKDAsMSxjb2w9Mixsd2Q9MikKICAgIGxlZ2VuZCgidG9wbGVmdCIscGFzdGUoImxhbWJkYT0iLGxhbWJkYSkpCn0KCiMjIyMjIyMjICBQV1MgcGFpcndpc2UgcnVuICMjIwpwd2Fwb3BzPC1jKCJQV1M5MSIgLCJQV1M5NiIsIlBXUzA3IiwgIlBXUzE3IikKY29tYjwtY29tYm4ocHdhcG9wcywgMikKY29tYjwtdChjb21iKQpjb21iMTwtY29tYltjKDEsNCw2LDMpLF0KCmZvciAoaSBpbiAxOm5yb3coY29tYjEpKXsKICAgIHBvcDE8LWNvbWIxW2ksMV0KICAgIHBvcDI8LWNvbWIxW2ksMl0KICAgIAogICAgI3M8LW5weUxvYWQocGFzdGUwKCJEYXRhL1BDQWFuZ3NkL3NlbGVjdGlvbi8iLHBvcDEsIi4iLHBvcDIsIl9zZWxlY3Rpb24uc2VsZWN0aW9uLm5weSIpKQogICAgczwtbnB5TG9hZChwYXN0ZTAoIi4uL0RhdGEvUENBYW5nc2Qvc2VsZWN0aW9uL3BydW5lZF8iLHBvcDEsIi4iLHBvcDIsIl9zZWxlY3Rpb24uc2VsZWN0aW9uLm5weSIpKQogICAgCiAgICAjaG93IG1hbnkgUEMgYXhlcyB3ZXJlIGV2YWx1YXRlZD8KICAgIG5jPC1uY29sKHMpCiAgICAKICAgICMjIG1ha2UgUVEgcGxvdCB0byBRQyB0aGUgdGVzdCBzdGF0aXN0aWNzCiAgICAjcXFjaGkocykKICAgIAogICAgIyBjb252ZXJ0IHRlc3Qgc3RhdGlzdGljIHRvIHAtdmFsdWUKICAgIGlmIChuYz09MSkgcHZhbDwtMS1wY2hpc3EocywxKQogICAgaWYgKG5jPjEpIHsKICAgICAgICBwJHB2YWwxPC1wdmFsWywxXQogICAgICAgIHAkcHZhbDI8LXB2YWxbLDJdCiAgICAgICAgcCRsb2M8LTE6bnJvdyhwKQogICAgICAgIHAkcHZhbDEubG9nPC0tbG9nMTAocCRwdmFsMSkKICAgICAgICBwJHB2YWwyLmxvZzwtLWxvZzEwKHAkcHZhbDIpCiAgICAgICAgfQogICAgCiAgICAjcmVhZCB0aGUgcG9zaXRpb24gaW5mbyAgICAKICAgIHA8LXJlYWQudGFibGUocGFzdGUwKCIuLi9EYXRhL1BDQWFuZ3NkL3NlbGVjdGlvbi9wcnVuZWRfIixwb3AxLCIuIixwb3AyLCJfc2VsZWN0aW9uLnNpdGVzIiksY29sQz1jKCJmYWN0b3IiLCJpbnRlZ2VyIiksc2VwPSI6IikKICAgIG5hbWVzKHApPC1jKCJjaHIiLCJwb3MiKQogICAgCiAgICAjIyBtYWtlIG1hbmhhdHRlbiBwbG90CiAgICBwJHB2YWw8LXB2YWwKICAgIHAkbG9jPC0xOm5yb3cocCkKICAgIHAkcHZhbC5sb2c8LS1sb2cxMChwJHB2YWwpCiAgICAKICAgIHdyaXRlLmNzdihwLCBwYXN0ZTAoIi4uL091dHB1dC9QQ0Evc2VsZWN0aW9uL0ZzdF9wdmFsX1BydW5lZCIsIHBvcDEsIi4iLHBvcDIsIi5jc3YiKSkKICAgIHdyaXRlLmNzdihwW3AkcHZhbC5sb2c+NCxdLCBwYXN0ZTAoIi4uL091dHB1dC9QQ0Evc2VsZWN0aW9uL1NlbGVjdGVkX2hpZ2hQX3NpdGVzX1BydW5lZCIsIHBvcDEsIi4iLHBvcDIsIi5jc3YiKSkKICAgIAogICAgY2g8LXVuaXF1ZShwJGNocikKICAgIAogICAgI2NvdW50IHRoZSBudW1iZXIgb2Ygc2l0ZXMgcGVyIGNocm9tb3NvbWVzCiAgICBwb3NzPC1kYXRhLmZyYW1lKGNocj1jaCkKICAgIGs9MQogICAgZm9yIChpIGluIDE6bnJvdyhwb3NzKSl7CiAgICAgICAgZGY8LXBbcCRjaHI9PXBvc3MkY2hyW2ldLF0KICAgICAgICBwb3NzJHN0YXJ0W2ldPC1rCiAgICAgICAgcG9zcyRlbmRbaV08LWsrbnJvdyhkZiktMQogICAgICAgIGs9aytucm93KGRmKQogICAgfQogICAgCiAgICBwb3NzJHg8LXBvc3Mkc3RhcnQrKHBvc3MkZW5kLXBvc3Mkc3RhcnQpLzIKICAgIAogICAgI2NvbG9yIHZlY3RvcnMKICAgIGNvbG9yczwtcmVwKGMoInN0ZWVsYmx1ZSIsImxpZ2h0Ymx1ZSIpLCB0aW1lcz0xMykKICAgIHAkY2hyPC1mYWN0b3IocCRjaHIsIGxldmVscz1wb3NzJGNocikKICAgIGdncGxvdChkYXRhPXAsIGFlcyh4PWxvYywgeT1wdmFsLmxvZywgY29sb3I9Y2hyKSkrCiAgICAgICAgZ2VvbV9wb2ludChzaXplPTAuMTUpKwogICAgICAgIHNjYWxlX2NvbG9yX21hbnVhbCh2YWx1ZXM9Y29sb3JzLCBndWlkZT0nbm9uZScpKwogICAgICAgIHNjYWxlX3hfY29udGludW91cyhuYW1lPSJDaHJvbW9zb21lIiwgYnJlYWtzPXBvc3MkeCwgbGFiZWxzPWdzdWIoImNociIsJycscG9zcyRjaHIpKSsKICAgICAgICB0aGVtZV9jbGFzc2ljKCkreWxhYigiLWxvZzEwKHAtdmFsdWUpIikrCiAgICAgICAgZ2d0aXRsZShwYXN0ZTAocG9wMSwiIHZzLiAiLHBvcDIpKQogICAgZ2dzYXZlKHBhc3RlMCgiLi4vT3V0cHV0L1BDQS9zZWxlY3Rpb24vcGNhbmdzZF9zZWxlY3Rpb25fIixwb3AxLCJfIixwb3AyLCJfcGxvdC5wbmciKSwgd2lkdGggPSA2LCBoZWlnaHQgPSAzLCBkcGkgPSAzMDApICAgIAp9CgpgYGAKIVtdKC4uL091dHB1dC9QQ0Evc2VsZWN0aW9uL3BjYW5nc2Rfc2VsZWN0aW9uX1BXUzkxX1BXUzk2X3Bsb3QucG5nKQohW10oLi4vT3V0cHV0L1BDQS9zZWxlY3Rpb24vcGNhbmdzZF9zZWxlY3Rpb25fUFdTOTZfUFdTMDdfcGxvdC5wbmcpICAKIVtdKC4uL091dHB1dC9QQ0Evc2VsZWN0aW9uL3BjYW5nc2Rfc2VsZWN0aW9uX1BXUzA3X1BXUzE3X3Bsb3QucG5nKQohW10oLi4vT3V0cHV0L1BDQS9zZWxlY3Rpb24vcGNhbmdzZF9zZWxlY3Rpb25fUFdTOTFfUFdTMTdfcGxvdC5wbmcpCgojIyBQQ0Fuc2dkIGZvciBhbGwgeWVhcnMgb2YgUFdTICAKYGBge3IgZXZhbD1GQUxTRSwgbWVzc2FnZT1GQUxTRSwgd2FybmluZz1GQUxTRX0KIyMgUFdTIHRvZ2V0aGVyCnM8LW5weUxvYWQoIi4uL0RhdGEvUENBYW5nc2Qvc2VsZWN0aW9uL1BXU19zZWxlY3Rpb24uc2VsZWN0aW9uLm5weSIpCgojIyBtYWtlIFFRIHBsb3QgdG8gUUMgdGhlIHRlc3Qgc3RhdGlzdGljcwojcXFjaGkocykKIyBjb252ZXJ0IHRlc3Qgc3RhdGlzdGljIHRvIHAtdmFsdWUKcHZhbDwtMS1wY2hpc3EocywxKQoKIyMgcmVhZCBwb3NpdGlvbnMgCnA8LXJlYWQudGFibGUoIi4uL0RhdGEvUENBYW5nc2Qvc2VsZWN0aW9uL1BXU19zZWxlY3Rpb24uc2l0ZXMiLGNvbEM9YygiZmFjdG9yIiwiaW50ZWdlciIpLHNlcD0iOiIpCm5hbWVzKHApPC1jKCJjaHIiLCJwb3MiKQpwJHB2YWwxPC1wdmFsWywxXQpwJGxvYzwtMTpucm93KHApCnAkcHZhbDEubG9nPC0tbG9nMTAocCRwdmFsMSkKCiNjb3VudCB0aGUgbnVtYmVyIG9mIHNpdGVzIHBlciBjaHJvbW9zb21lcwpjaDwtdW5pcXVlKHAkY2hyKQpwb3NzPC1kYXRhLmZyYW1lKGNocj1jaCkKaz0xCmZvciAoaSBpbiAxOm5yb3cocG9zcykpewogICAgZGY8LXBbcCRjaHI9PXBvc3MkY2hyW2ldLF0KICAgIHBvc3Mkc3RhcnRbaV08LWsKICAgIHBvc3MkZW5kW2ldPC1rK25yb3coZGYpLTEKICAgIGs9aytucm93KGRmKQp9Cgpwb3NzJHg8LXBvc3Mkc3RhcnQrKHBvc3MkZW5kLXBvc3Mkc3RhcnQpLzIKCmNvbG9yczwtcmVwKGMoInN0ZWVsYmx1ZSIsImxpZ2h0Ymx1ZSIpLCB0aW1lcz0xMykKcCRjaHI8LWZhY3RvcihwJGNociwgbGV2ZWxzPXBvc3MkY2hyKQoKZ2dwbG90KGRhdGE9cCwgYWVzKHg9bG9jLCB5PXB2YWwxLmxvZywgY29sb3I9Y2hyKSkrCiAgICBnZW9tX3BvaW50KHNpemU9MC4xKSsKICAgIHNjYWxlX2NvbG9yX21hbnVhbCh2YWx1ZXM9Y29sb3JzLCBndWlkZT0nbm9uZScpKwogICAgc2NhbGVfeF9jb250aW51b3VzKG5hbWU9IkNocm9tb3NvbWUgcG9zaXRpb24iLCBicmVha3M9cG9zcyR4LCBsYWJlbHM9Z3N1YigiY2hyIiwnJyxwb3NzJGNocikpKwogICAgdGhlbWVfY2xhc3NpYygpK3lsYWIoIi1sb2cxMChwLXZhbHVlKSIpKwogICAgZ2d0aXRsZSgiUFdTIikKZ2dzYXZlKCIuLi9PdXRwdXQvUENBL3NlbGVjdGlvbi9QV1NfcGNhbmdzZF9zZWxlY3Rpb24ucG5nIiwgd2lkdGggPSAxMiwgaGVpZ2h0ID0gNiwgZHBpPTMwMCkgICAgCmBgYAohW10oLi4vT3V0cHV0L1BDQS9zZWxlY3Rpb24vUFdTX3BjYW5nc2Rfc2VsZWN0aW9uLnBuZykKCiMjIEhpZ2ggUC12YWx1ZSBsb2NpIGZyb20gZWFjaCBwb3AgcGFpcgoKYGBge3IgZWNobz1UUlVFLCBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQojZm9yIGVhY2ggcG9wdWxhdGlvbiBwYWlyCnB3YXBvcHM8LWMoIlBXUzkxIiAsIlBXUzk2IiwiUFdTMDciLCAiUFdTMTciKQpjb21iPC1jb21ibihwd2Fwb3BzLCAyKQpjb21iPC10KGNvbWIpCmNvbWIxPC1jb21iW2MoMSw0LDYsMyksXQoKI2NyZWF0ZSB3aW5kb3dzIHRvIGJlIGFzc2lnbmVkIGZvciBlYWNoIHNpdGUKY2hyIDwtIHJlYWQudGFibGUoJy4uL0RhdGEvbmV3X3ZjZi9jaHJfc2l6ZXMuYmVkJykKY2hyPC1jaHJbLC0yXQpjb2xuYW1lcyhjaHIpPC1jKCJjaHIiLCAibGVuIikKI2Nocm1heCRzdGFydDwtYygwLGN1bXN1bShjaHJtYXgkbGVuKVsxOihucm93KGNocm1heCktMSldKQojc2V0a2V5KGNocm1heCwgY2hyKQoKd2luc3ogPC0gNTAwMDAgIyB3aW5kb3cgc2l6ZQojIHVzZSBzZXEgdG8gZmluZCB0aGUgc3RhcnQgcG9pbnRzIG9mIGVhY2ggd2luZG93CndpbmRvd3M8LWRhdGEuZnJhbWUoKQpmb3IgKGkgaW4gMTpucm93KGNocikpewogICAgY2g8LWNociRsZW5bY2hyJGNocj09cGFzdGUwKCJjaHIiLGkpXQogICAgd2luZG93X3N0YXJ0IDwtIHNlcShmcm9tID0gMSwgdG8gPSBjaC01MDAwMCwgYnkgPSB3aW5zeikKICAgIHdpbmRvd19zdG9wIDwtIHdpbmRvd19zdGFydCArICh3aW5zei0xKQogICAgd2luIDwtIGRhdGEuZnJhbWUoc3RhcnQgPSB3aW5kb3dfc3RhcnQsIHN0b3AgPSB3aW5kb3dfc3RvcCwgCiAgICAgICAgICAgICAgICAgICAgICBtaWQgPSB3aW5kb3dfc3RhcnQgKyAoKHdpbmRvd19zdG9wLXdpbmRvd19zdGFydCkvMi0wLjUpKQogICAgd2luJGNocjwtcGFzdGUwKCJjaHIiLGkpCiAgICB3aW4kY2g8LWkKICAgIHdpbmRvd3M8LXJiaW5kKHdpbmRvd3MsIHdpbikKfQpzZXREVCh3aW5kb3dzKQpmc3RkYXQ8LWxpc3QoKQpGc3Q8LWRhdGEuZnJhbWUoKQpmb3IgKGkgaW4gMTogbnJvdyhjb21iMSkpewogICAgcG9wMTwtY29tYjFbaSwxXQogICAgcG9wMjwtY29tYjFbaSwyXQogICAgZGY8LWZyZWFkKHBhc3RlMCgiLi4vT3V0cHV0L1BDQS9zZWxlY3Rpb24vU2VsZWN0ZWRfaGlnaFBfc2l0ZXNfUHJ1bmVkIiwgcG9wMSwiLiIscG9wMiwiLmNzdiIpKQogICAgZGY8LWRmWywtMV0KICAgICNhc3NpZ24gdGhlIHdpbmRvdyBpbmZvCiAgICBkZjwtd2luZG93c1tkZiwgLihjaHIsIHBvcywgcHZhbCwgcHZhbC5sb2csIHN0YXJ0LCBzdG9wLCBtaWQpLCBvbj0uKGNocj1jaHIsIHN0YXJ0IDw9cG9zLCBzdG9wPj1wb3MpXQogICAgI2NyZWF0ZSBhIHVuaXF1ZSBpZCBmb3IgdGhlIHdpbmRvdwogICAgZGYkd2luX2lkPC1wYXN0ZTAoZGYkY2hyLCJfIixkZiRtaWQpCiAgICBkZiRwb3BzPC1wYXN0ZTAocG9wMSwiLiIscG9wMikKICAgIGZzdGRhdFtbaV1dPC1kZgogICAgbmFtZXMoZnN0ZGF0KVtpXTwtcGFzdGUwKHBvcDEsIi4iLHBvcDIpCiAgICBGc3Q8LXJiaW5kKEZzdCwgZGYpCn0KCiNGaW5kIG92ZXJsYXBwaW5nIHBvc2l0aW9ucyBhbW9uZyBwYWlycz8KbGFwcGx5KGZzdGRhdCwgIltbIiwgIndpbl9pZCIpCgpSZWR1Y2UoaW50ZXJzZWN0LCBsaXN0KGZzdGRhdFtbMV1dJHdpbl9pZCAsIGZzdGRhdFtbMl1dJHdpbl9pZCkpClJlZHVjZShpbnRlcnNlY3QsIGxpc3QoZnN0ZGF0W1syXV0kd2luX2lkICwgZnN0ZGF0W1szXV0kd2luX2lkKSkKUmVkdWNlKGludGVyc2VjdCwgbGlzdChmc3RkYXRbWzNdXSR3aW5faWQgLCBmc3RkYXRbWzRdXSR3aW5faWQpKQoKRnN0JGNocjwtZmFjdG9yKEZzdCRjaHIsIGxldmVscz1wYXN0ZTAoImNociIsMToyNikpCmdncGxvdChGc3QsIGFlcyh4PXBvcywgeT1wdmFsLmxvZywgY29sb3I9cG9wcykpKwogICAgZ2VvbV9wb2ludCgpKwogICAgZmFjZXRfd3JhcCh+Y2hyKSsKICAgdGhlbWVfYncoKSt4bGFiKCJHZW5vbWUgcG9zaXRpb24iKSsKICAgIHlsYWIoIi1sb2cxMChQLXZhbCkiKSt0aGVtZShsZWdlbmQudGl0bGUgPSBlbGVtZW50X2JsYW5rKCksIGF4aXMudGV4dC54ID0gZWxlbWVudF9ibGFuaygpKQpnZ3NhdmUoIi4uL091dHB1dC9Gc3QvcGNhbmdzZF9zY2FuX2hpZ2hQc2l0ZXMucG5nIiwgd2lkdGg9OSwgaGVpZ2h0PTYsIGRwaT0zMDApCmBgYAohW10oLi4vT3V0cHV0L0ZzdC9wY2FuZ3NkX3NjYW5faGlnaFBzaXRlcy5wbmcpCgoKCjxicj4KPGJyPgoKCg==